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Abstract. We present new 3D magnetohydrodynamic (MHD) simulations of a supern ova- 
driven, stratified interstellar medium. These simulations were run using the Waagan 



et al. ( |2011| ) positivity preserving scheme for ideal MHD implemented in the Flash 
code. The scheme is stable even for the Mach numbers approaching 100 found in this 
problem. We have previously shown that the density distribution arising from hydro- 
dynamical versions of these simulations creates low-density pathways through which 
Lyman continuum photons can travel to heights \z\ > 1 kpc. This naturally produces the 
warm ionized medium through photoionization due primarily to O stars near the plane. 
However, our earlier models reproduce the peak but not the width of the observed emis- 
sion measure distribution. Here, we examine whether inclusion of magnetic fields and 
a greater vertical extent to the simulation domain produce a gas distribution that bet- 
ter matches the observations. We further study the change of magnetic energy over 
time in our models, showing that it appears to reach a steady state after a few hun- 
dred megayears, presumably supported by a turbulent dynamo driven by the supernova 
explosions. 
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1. Observational Motivation 



Supernova explosions provide one of the major sources of energy for observed tur- 
bulent motions in the interstellar gas. It can be argued that in star forming regions 
of galaxies their influence predominates ( |Mac Low & Klessen||2004[ |Tamburro et al. 
|2009| ). To constrain their influence on the interstellar gas, we have compared observa- 
tions of the upper atmosphere and gaseous halo of the Galaxy to numerical simulations 
of supernova-driven turbulence. 

The base model that we used, descr ibed by|Joung & Mac Low| ( |2006] ) and |Joung 
et al.|(2009), and inspired by the work of |Avillez] ^ 2000 ), used version 2.4 of the Flash 



code ( Fryxell et al.||2000) to model a square section through a stratified galactic disk 
centered on the midplane, with dimensions 0.5 x 0.5 x 10 kpc and periodic boundary 
conditions in the horizontal direction. The vertical boundary conditions were taken 
to be constant pressure to allow outflow. The adaptive mesh refinement capability of 
Flash was used to fix the maximum resolution within 200 pc of the midplane to 1.95 pc, 
with reduced maximum resolution at higher altitudes for computation efficiency. Super- 
nova explosions were placed at both random and correlated locations to reproduce their 
distribution in the Galaxy. Radiative cooling was implemented with an equilibrium 
ionization, solely temperature dependent cooling function, and a diffuse heating term 
characteristic of photoelectric heating and stratified with altitude was also included. 

|Henley et al.| ( |2010| ) measured X-ray emission from the hot, gaseous halo of the 
Galaxy by using shadowing against relatively nearby clouds to remove contributions 
from local sources including solar wind charge exchange and emission from the Local 
Bubble. They then compared these observations to simulations of emission from our 
models. They found good agreement between the simulated and observed emission 
measure of hot gas looking out of the Galaxy. However, the effective temperature of 
the emission in the model systematically climbed with time. 

| Wood et aL] ( |2010| ) compared the distribution of Ha emission measures measured 
by the Wisconsin Ha Mapper (WHAM Haffne r~et al.|20 03) along lines of sight out of 
the plane of the Galaxy to the distribution of simulated Ha emission from the model (see 
Fig. [TfaJ). Figure [TfZ? ) shows the comparison. They found that the emission measure 
distribution of the models was centered around the same value as the observations, but 
that it was broader, with more lines of sight having very high or low emission measures 
seen in the models. This suggested that some additional homogenizing effect might be 
acting beyond the physics included in the models. Magnetic fields could have such an 
effect, so we decided to investigate how their inclusion would change our results. We 
here report on our implementation of fields in our models and some first results, though 
we only reach preliminary conclusions. 



2. Simulation Techniques 

Our problem requires the computation of flows with Alfven and Mach numbers greatly 
exceeding unity in plasmas with a range of plasma J3 (ratio of thermal to magnetic 
pressure) extending from well under to well over unity. Most linearized Riemann 
solvers have substantial stability problems under these conditions. Therefore we used 
the positivity-preserving generalization by |Bouchut et al.| ( |2007| ) of an HLL (Harten 



et al.||1983) approximate Riemann solver with a base MUSCL-Hancock scheme. This 



solver satisfies discrete entropy inequalities and imposes positivity of density and in- 
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Figure 1. (a) Edge-on visualizations of em ission measures ( J n^ds) from a sim- 
ulation grid described by |Joung et al. ( 2009 ). Adapted from Fig. 3 of | Wood et al. 
(2010). (b) Emission measure distributions for the Galactic DIG from the WHAM 
survey and from our photoionization simulation. Fraction (dex) -1 is the fraction of 
points in each bin divided by the logarithmic bin interval. The black histogram is of 
J n 2 e ds in the simulation, viewed face-on. The red histogram is of EM sin \b\ from 
the WHAM survey with classical H n regions and sightlines with \b\ < 10° removed, 
from [Hill et alT] ([2008 ). Lognormal fits to each distribution are also shown. For 
the emission measures from the simulations we have rem oved contributions from 
midplane regions with \z\ < 150 pc. Adapted from Fig. 6 of Woo det al.| ( |2010] ). 



ternal energy. The three-wave HLL3R version of |Bouchut et al.| ( [2010| ) is used here, in 
an implementation described by Waagan (2009]). This solver vastly improves stability, 
and thus efficiency, without sacrificing accuracy. 

Waa gan et al.| ( |2011| ) demonstrates that HLL3R performs as well as other HLL 



schemes on standard test problems such as shock tubes. Figure [2] shows a comparison 
of solutions of an Orszag-Tang vortex by two of the standard Flash solvers, Roe and 
HLLE to the HLL3R solution. The accuracy compares to the Roe solver, while the 
slightly more diffusive HLLE solver fails to capture the central structure. However, 
the comparison of the performance of the Roe solver to HLL3R for driven Mach 2 
turbulence given in Table 3 of Wa agan et al.| ( |201 \\ shows that the Roe solver requires 
a Courant-Friedrich-Lewy number below 0.2 for stability, while HLL3R remains stable 
at 0.8. Even when the Roe solver doesn't crash outright, it regularly generates hot 
zones that force the Courant timestep to low values. As a result, the total CPU time 
used for this problem is a factor of 9 lower for HLL3R. Increasing the Mach number to 
10 in J3 = 1 plasma resulted in a problem that simply could not be stably computed by 
the Roe or HLLE solvers, but was solved by HLL3R. This stability was demonstrated 



under even more extreme conditions in the driven MHD turbulence models by Brunt 
et al. (2010) which reached Mach numbers as high as 20 with J3 = 1.25 x 10" 3 . . 



Our magnetized, supernova-driven models were initialized with constant /? fields 
having midplane strengths of either 6.5 or 13 fiG. We found during initial testing that 
the grid we used in Joung & Mac Low| ( |2006[ ) and |Joung et al.| ([2009) was insufficiently 
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Figure 2. Orszag-Tang test. Pressure at time t = 1 with resolution h = 256 1 . 
Top to bottom left to right: FLASH-Roe, HLL3R, FLASH-HLLE. Figure 4.3 from 
Waaganet aL]([2QTT1) 



large to capture the dynamics of the halo, though it sufficed to model the midplane 
gas that we focused on in those papers. Instead, we used 1 x 1 x 40 kpc grids for our 
standard models. Also, instead of initializing an entirely isothermal grid at 10 4 K, we 
set the halo temperature to 10 6 K at altitudes \z\ > 1 kpc to avoid very strong shocks at 
high altitude during the initial expansion of supernova blast waves. 



3. Preliminary Results 



3.1. Vertical Structure 

The magnetic field indeed leads to somewhat smoother flows, with more filamentary 
structure seen above a kiloparsec. However, the biggest surprise came from examin- 
dation of the vertical structure of the halo. The ever increasing halo temperature seen 



by |Henley et al.| ( |2010| ) was shown to be a numerical artifact. The constant-pressure 
vertical boundary conditions used by Joun g & Mac Low| ( |2006| ) and |Joung et al.| ([2009) 
did not have a velocity limiter to prevent inflow. After a hot initial transient flow from 
the first supernova explosions in the halo crossed the boundary, the pressure and tem- 
perature were set high. The gas then started to fall inward lacking further support, and 
the boundary then fed high pressure and temperature continuously on to the grid. Rais- 
ing the boundary and setting it to allow only outflow resulted in much more reasonable 
behavior, with an appropriately thick disk not confined by high pressure gas, and a 
moderate temperature halo, as shown in Figure [3] In future work we will examine in 
detail whether the effective X-ray temperature, Ha emission measure, and H i 21 cm 
distributions now reproduce the observations. 

However, we then saw that even in steady state, the halo undergoes vertical oscil- 
lation on either side of the midplane, uncorrelated between the upper and lower sides. 
Exactly this behavior was predicted by | Walters & Cox| ( |2001| ) using a one-dimensional 
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Figure 3. Slices of temperature over time for a model with initial midplane mag- 
netic field of 6.5 jiG. The time of each slice in megayears is listed above each image. 
Note that the aspect ratio is not unity; the vertical scale of these images is contracted. 



Adapted from Figure 2 of |Hill et al. (2011 



model and proposed as an explanation for the remarkable observation that gas towards 
the northern Galactic pole seems to be primarily infalling while that towards the south- 
ern Galactic pole seems to primarily flow outwards. Although our models reproduce 
only a small portion of a galactic disk, we may take from them the qualitative statement 
that the upper atmospheres and halos of galaxies are unlikely to be quiescent, static 
places, but rather are likely to be in continuous hypersonic motion, with observable 
consequences. 



3.2. Dynamos 



The modeling of an idealized turbulent dynamo using HLL3R was already demon- 
strated by Waaga net al.| ( |20ir] ). This work was extended to consider dynamos in ac- 
cretion flows with application to the formation of the first stars by |Sur et al.] ( |2010| ). 
Federrath et al.| ( |2011b| ) demonstrated that dynamo action in accretion only occurs for 
resolutions of at least 30 cells per Jeans length, and does not show convergence out 
to at least 128 cells per Jeans length, instead showing continuously stronger growth of 



rms field with higher resolution and higher effective Reynolds number (Federrath et al. 
201 la| ). |Balsara et al.| ( |2004| ) had already shown that supernova-driven turbulence in a 
periodic box acts as a turbulent dynamo. |Hill et al.] ( |2011| ) shows that this likely holds 
in a stratified medium as well, with magnetic field energies initially varying dependent 
on their initial values but then appearing to converge towards values characteristic of 
field strengths somewhat under 5 /iG. However runs of 400 Myr still had not clearly 
moved into a steady-state regime, so this conclusion remains preliminary. 
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4. Conclusions 



Supernova-driven turbulence appears increasingly likely to be able to explain observa- 
tions in Ha and X-ray of material above the plane of the Galaxy, along with many other 
observations. Hypersonic vertical motions of large regions of galactic halos above star- 
forming disks appear to be driven by transient variations in energy input in the much 



denser disk midplanes of the galaxies, as predicted by [Walters & Cox| ( |2001| ). These 
vertical motions appear to explain at least some of the discrepancies between models 
and observations identified by He nley et"aL] ( |2010| ) and |Wood et al.| ( |20T0| ). Our mag- 
netized models demonstrated that magnetic fields can be modeled stably even in highly 
compressible, highly magnetized flows using the positivity preserving HLL3R solver 
described by Waagan (2009). They further suggested that a turbulent dynamo can be 
driven in a stratified disk, as well as in a periodic box. 
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